﻿using System;
using System.IO;
using VCPPDiffusion2D;

namespace NumericalSolution {
  class VCPPDiffusion2DMain {
    static int numberOfNodesInX = 101;
    static int numberOfNodesInY = 101;
    static double deltaT = 0.01;
    static int endCount = 200;
    static bool outswitch = false;

    public static void Main(String[] args) {

      VCPPDiffusion2DADI sim = new VCPPDiffusion2DADI();
      sim.NumberOfNodesInX = numberOfNodesInX;
      sim.NumberOfNodesInY = numberOfNodesInY;
      sim.DeltaX = 1.0 / (numberOfNodesInX - 1);
      sim.DeltaY = 1.0 / (numberOfNodesInY - 1);
      sim.DeltaT = deltaT;
      sim.Alpha = 0.005;

      double[] u = new double[numberOfNodesInX * numberOfNodesInY];
      PrepareInitialValues(sim, u);

      for (int i = 0; i <= endCount; i++) {
        if (outswitch == true) {
          double timeInMiliseconds = i * sim.DeltaT * 1000.0;
          Print(sim, u, "c:/nsworkspace/csvcppadi2d_"
            + timeInMiliseconds.ToString("00000") + ".txt");
        }
        sim.NextStep(u);
      }
    }

    public static void PrepareInitialValues(VCPPDiffusion2DADI sim, 
                                            double[] u) {

      for (int j = 0; j < sim.NumberOfNodesInY; j++) {
        for (int i = 0; i < sim.NumberOfNodesInX; i++) {
          double x = i * sim.DeltaX;
          double y = j * sim.DeltaY;

          u[j * sim.NumberOfNodesInX + i] = 
            Math.Exp(-0.5 * Math.Pow((x - 0.5) / 0.1, 2.0)
                     - 0.5 * Math.Pow((y - 0.5) / 0.1, 2.0));

          if (i == 0 || i == numberOfNodesInX - 1 ||
              j == 0 || j == numberOfNodesInY - 1) {
            u[j * sim.NumberOfNodesInX + i] = 0.0;
          }
        }
      }
    }

    public static void Print(VCPPDiffusion2DADI sim, 
                             double[] u, String outfile) {

      StreamWriter writer = File.CreateText(outfile);
      for (int i = 0; i < sim.NumberOfNodesInX; i++) {
        for (int j = 0; j < sim.NumberOfNodesInY; j++) {
          double x = i * sim.DeltaX;
          double y = j * sim.DeltaY;
          writer.WriteLine(x + " " + y + " " 
            + u[j * sim.NumberOfNodesInX + i]);
        }
        writer.WriteLine();
      }
      writer.Close();
    }
  }

}
